library(tidyverse)
library(gifski)
library(ggraph)
library(here)
library(igraph)
source(here("modelFunction_rewiring.R"))
First, run the model.
# Define parameters
N = 100
edge.prob <- 0.04
burn.in = 20
burn.out = 5
pm = 0.3
ps = 0.1
pa = 0.2
add00 = c(0.5, 10)
lose01 = 0.1
add10 = 0.05
lose11 = c(0.5, 0.5)
histMultiplier = 1.2
doRemoval = TRUE
modelGraphs <- runModel(N = N, # Nodes in the network
edge.prob = edge.prob,
burn.in = burn.in,
burn.out = burn.out,
pm = pm,
ps = ps,
pa = pa,
add00 = add00,
lose01 = lose01,
add10 = add10,
lose11 = lose11,
histMultiplier = histMultiplier,
doRemoval = doRemoval) %>%
lapply(., function(x){
igraph::graph_from_adjacency_matrix(x, mode = "undirected")
})
Compute degree for each of the model networks.
degrees <- lapply(modelGraphs, function(x){igraph::degree(x)})
degreesDF <- do.call(cbind, degrees) %>% as.data.frame() %>% setNames(., 1:length(degrees)) %>% mutate(individual = 1:nrow(.)) %>%
pivot_longer(cols = -individual, names_to = "slice", values_to = "degree") %>%
mutate(slice = as.numeric(slice)) %>%
mutate(beforeAfter = ifelse(slice <= burn.in, "before", "after"))
## Warning in (function (..., deparse.level = 1) : number of rows of result is not
## a multiple of vector length (arg 21)
Okay, so we get some degree distributions, but do we get emergent individual variation in degree that’s somewhat consistent throughout time?
# Plot everyone's degree over time
degreesDF %>%
ggplot(aes(x = slice, y = degree, col = individual, group = individual))+
geom_line()+
theme_minimal()+
geom_vline(xintercept = burn.in+1, col = "red", size = 2)
# Let's just look at a few individuals, because this plot is too hard to read
indivs <- sample(1:max(degreesDF$individual), 5)
degreesDF %>%
filter(individual %in% indivs) %>%
ggplot(aes(x = slice, y = degree, col = factor(individual), group = individual))+
geom_line()+
theme_minimal()+
geom_vline(xintercept = burn.in+1, col = "red", size = 2)
# This is still a little hard to see... let's look at just the "before" rows.
degreesDF %>%
filter(individual %in% indivs, beforeAfter == "before") %>%
ggplot(aes(x = slice, y = degree, col = factor(individual), group = individual))+
geom_line()+
theme_minimal()
# Okay, I'm seeing some differentiation, but no clear patterns. Because of how the model is set up, the rank order tends to be maintained for a few consecutive time steps, but not over the whole span of the baseline model. This makes sense, since individuals don't have any a priori reason to have higher or lower degrees than other individuals.
Now, going to compute three network-level measures before and after the removal of an individual.
How does edge density behave over time?
densities <- lapply(modelGraphs, function(x){
edge_density(x)
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "density") %>%
mutate(slice = 1:length(modelGraphs))
densities %>%
ggplot(aes(x = slice, y = density))+
geom_line()+
geom_smooth()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
First delta: immediately before removal to after removal/pre-rewiring
delta1_density <- igraph::edge_density(modelGraphs[[burn.in+1]])-igraph::edge_density(modelGraphs[[burn.in]])
Second delta: after removal/pre-rewiring to post-rewiring
delta2_density <- igraph::edge_density(modelGraphs[[burn.in+2]])-igraph::edge_density(modelGraphs[[burn.in+1]])
How does connectivity behave over time?
connectivities <- lapply(modelGraphs, function(x){
mean_distance(x)
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "connectivity") %>%
mutate(slice = 1:length(modelGraphs))
connectivities %>%
ggplot(aes(x = slice, y = connectivity))+
geom_line()+
geom_smooth()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 row(s) containing missing values (geom_path).
First delta: immediately before removal to after removal/pre-rewiring
delta1_meanDistance <- igraph::mean_distance(modelGraphs[[burn.in+1]])-igraph::mean_distance(modelGraphs[[burn.in]])
Second delta: after removal/pre-rewiring to post-rewiring
delta2_meanDistance <- igraph::mean_distance(modelGraphs[[burn.in+2]])-igraph::mean_distance(modelGraphs[[burn.in+1]])
How does modularity behave over time?
modularities <- lapply(modelGraphs, function(x){
cluster_fast_greedy(x) %>%
modularity()
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "modularity") %>%
mutate(slice = 1:length(modelGraphs))
modularities %>%
ggplot(aes(x = slice, y = modularity))+
geom_line()+
geom_smooth()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
First delta: immediately before removal to after removal/pre-rewiring
delta1_modularity <- modularity(cluster_fast_greedy(modelGraphs[[burn.in+1]]))-modularity(cluster_fast_greedy(modelGraphs[[burn.in]]))
Second delta: after removal/pre-rewiring to post-rewiring
delta2_modularity <- modularity(cluster_fast_greedy(modelGraphs[[burn.in+2]]))-modularity(cluster_fast_greedy(modelGraphs[[burn.in+1]]))
I am not at all convinced that these momentary deltas are what I should be using. Seems like we need to at least consider two time steps out, since the model relies on this. Looking at the graphs above shows me how different the longer-term trends are from the short-term fluctuations.
Let’s run the model 100 times and see if trends emerge in how to measure.
nRuns <- 100
modelRuns <- vector(mode = "list", length = nRuns)
for(i in 1:nRuns){
modelRuns[[i]] <- runModel(N = N, # Nodes in the network
edge.prob = edge.prob,
burn.in = burn.in,
burn.out = burn.out,
pm = pm,
ps = ps,
pa = pa,
add00 = add00,
lose01 = lose01,
add10 = add10,
lose11 = lose11,
histMultiplier = histMultiplier,
doRemoval = doRemoval) %>%
lapply(., function(x){
igraph::graph_from_adjacency_matrix(x, mode = "undirected")
})
}
densList <- lapply(modelRuns, function(x){
lapply(x, function(y){
edge_density(y)
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "density") %>%
mutate(slice = 1:nrow(.))
})
densDF <- data.table::rbindlist(densList, idcol = "run")
connList <- lapply(modelRuns, function(x){
lapply(x, function(y){
mean_distance(y)
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "connectivity") %>%
mutate(slice = 1:nrow(.))
})
connDF <- data.table::rbindlist(connList, idcol = "run")
moduList <- lapply(modelRuns, function(x){
lapply(x, function(y){
cluster_fast_greedy(y) %>%
modularity()
}) %>%
unlist() %>%
as.data.frame() %>%
setNames(., "modularity") %>%
mutate(slice = 1:nrow(.))
})
moduDF <- data.table::rbindlist(moduList, idcol = "run")
Make some plots and try to detect any changepoint trends
densDF %>%
ggplot(aes(x = slice, y = density))+
geom_line(aes(group = run, col = run), size = 0.1)+
geom_smooth()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
scale_color_viridis_c()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
connDF %>%
ggplot(aes(x = slice, y = connectivity))+
geom_line(aes(group = run, col = run), size = 0.1)+
geom_smooth()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
scale_color_viridis_c()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 100 rows containing non-finite values (stat_smooth).
## Warning: Removed 100 row(s) containing missing values (geom_path).
moduDF %>%
ggplot(aes(x = slice, y = modularity))+
geom_line(aes(group = run, col = run), size = 0.1)+
geom_smooth()+
theme_minimal()+
geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
scale_color_viridis_c()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'